A light introduction to runoutSIM

Jason Goetz

2025-08-19

Introduction

The runoutSIM package is developed for regional-scale runout simulations (spatial distribution and velocity) using random walks. This vignette provides an introduction to applying runoutSIM for regional debris flow simulation using a subcatchment of the upper Maipo river basin, Chile as an example. In this vignette, we will:

Loading packages and reading spatial data

Spatial handling in the runoutSIM package is supported by the sf (vector data) and terra (gridded/raster data) packages. In this example, we load/read a 12.5 m spatial resolution DEM, debris flow runout polygons, and corresponding source points mapped for each polygon. We also create a hillshade model to help with visualization of the results.

Tip: Currently, this package is designed to work with data in a UTM coordinate reference system (CRS). You will need to project your own data to a local UTM CRS for it to work.

# load packages
library(runoutSim)
library(terra)
#> Warning: package 'terra' was built under R version 4.4.3
#> terra 1.8.54
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

# Load digital elevation model (DEM)
dem <- rast("C:/GitProjects/runoutSIM/Dev/Data/elev_fillsinks_WangLiu.tif")

# Compute hillshade for visualization 
slope <- terrain(dem, "slope", unit="radians")
aspect <- terrain(dem, "aspect", unit="radians")
hill <- shade(slope, aspect, 40, 270)

# Load debris flow runout source points and polygons
source_points <- st_read("C:/GitProjects/runoutSIM/Dev/Data/debris_flow_source_points.shp")
#> Reading layer `debris_flow_source_points' from data source 
#>   `C:\GitProjects\runoutSim\Dev\Data\debris_flow_source_points.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 73 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 389175.6 ymin: 6293926 xmax: 398788.6 ymax: 6327439
#> Projected CRS: WGS 84 / UTM zone 19S
runout_polygons <- st_read("C:/GitProjects/runoutSIM/Dev/Data/debris_flow_runout_polygons.shp")
#> Reading layer `debris_flow_runout_polygons' from data source 
#>   `C:\GitProjects\runoutSim\Dev\Data\debris_flow_runout_polygons.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 73 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 389139.1 ymin: 6293864 xmax: 398852.9 ymax: 6327455
#> Projected CRS: WGS 84 / UTM zone 19S

# Load basin boundary
bnd_catchment <- st_read("C:/GitProjects/runoutSIM/Dev/Data/basin_rio_olivares.shp")
#> Reading layer `basin_rio_olivares' from data source 
#>   `C:\GitProjects\runoutSim\Dev\Data\basin_rio_olivares.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 15 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 386042.8 ymin: 6293348 xmax: 404405.3 ymax: 6331286
#> Projected CRS: WGS 84 / UTM zone 19S

Creating interactive maps

Using the leafmap() function, we can quickly create an interactive leaflet map for viewing the data. Similar to dplyr, we use the %>% operator to pipe together multiple leafmap() calls to create a combined leaflet map with multiple layers.

Tip: leafmap() has the option to interactively query raster values with a mouse hover. This option is turned on by default, but can result in a considerable increase in the file size of the leaflet widget (e.g. when exported for sharing). It is recommended to have add_image_query = FALSE when saving / exporting the leaflet map.

map <- leafmap(bnd_catchment, color = '#f7f9f9', fill_color = '#FF000000', weight = 4) %>%
    leafmap(runout_polygons) %>%
    leafmap(source_points, color = '#e74c3c') %>%
    leafmap(hill, palette = grey(0:100/100), opacity = 1, add_legend = FALSE, 
            add_image_query = FALSE) %>%
    leafmap(dem, palette = viridis::mako(10), opacity = 0.6, add_image_query = FALSE)
map

Simulate runout from a source cell

Runout paths are simulated from individual source cells, defined as x, y coordinates. For example, if we are simulating runout from a single source cell defined by a source point (vector data), we use sf::st_coordinates() to extract the x, y values.

# Select a single debris flow and source point for the example
runout_polygon <- runout_polygons[31,]

# Get corresponding source point
source_point  <- st_filter(st_as_sf(source_points), st_as_sf(runout_polygon))

# Get coordinates of source point
source_crds <- st_coordinates(source_point)
print(source_crds)
#>             X       Y
#> [1,] 395698.1 6307421

Now that we have the coordinates, we can run the random walk runout simulations using runoutSim(). Where:

sim_paths = runoutSim(dem = dem, 
                      xy = source_crds, 
                      mu = 0.08, 
                      md = 140, 
                      slp_thresh = 35, 
                      exp_div = 2.5, 
                      per_fct = 1.95, 
                      walks = 1000)

# plot structure of sim_paths
str(sim_paths)
#> List of 4
#>  $ start_cell    : num 2803625
#>  $ cell_trav_freq: 'table' int [1:2012(1d)] 5 2 7 17 12 7 4 1 1 1 ...
#>   ..- attr(*, "dimnames")=List of 1
#>   .. ..$ all_indices: chr [1:2012] "2774166" "2774167" "2775633" "2775634" ...
#>  $ cell_max_vel  : num [1:2012] 15.8 16.2 12 14.3 16.2 ...
#>  $ prob_connect  : NULL

runoutSim() returns a list containing:

These values are stored as cell numbers (e.g., $all_indices). Although not calculated here, prob_connect provides the probability the source cell connected to a feature of interest.

Converting results to a raster

To visualize the results, we convert the output from runoutSim() to a raster.

## Convert paths to raster with cell transition frequencies
paths_raster <- walksToRaster(sim_paths, dem, method = "freq")

# Plot results
paths_raster <- crop(paths_raster, ext(runout_polygon)+500)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Runout traverse frequency")
plot(paths_raster, legend = T, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)

In walksToRaster(), we use method = "freq" to return traverse frequencies. We can also return traverse (ECDF) probabilities by using method = "cdf_prob".

We use velocityToRaster() to get a raster of velocity values (m/s).

## Convert paths to raster with max. velocities
vel_raster <- velocityToRaster(sim_paths, dem)

# Plot results
vel_raster <- crop(vel_raster, ext(runout_polygon)+500)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Runout velocity")
plot(vel_raster, col = map.pal('plasma'), legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)

Computing connectivity probability

The runoutSim() function also allows us to calculate the connectivity probability of a source cell to feature of interest - i.e. the proportion of walks that intersect a feature from a given source cell.

In this example, we will determine the connectivity probability of a simulated debris flow to a river channel. Using makeConnFeature() function, we will create a matrix indicating with a value of 1 which cells contain the river channel.

# Load river channel polygon (vector)
river_channel <- st_read("C:/GitProjects/runoutSIM/Dev/Data/river_channel.shp")
#> Reading layer `river_channel' from data source 
#>   `C:\GitProjects\runoutSim\Dev\Data\river_channel.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 393504.8 ymin: 6293330 xmax: 396806.5 ymax: 6328149
#> Projected CRS: WGS 84 / UTM zone 19S

# Create a connectivity feature for runoutSim that matches the input DEM
feature_mask <- makeConnFeature(river_channel, dem)

Now that we have a connectivity feature, we can re-run runoutSIM() with some additional parameter options: source_connect = TRUE allowing to calculate a connectivity probability and connect_feature = feature_mask.

sim_paths = runoutSim(dem = dem, 
                      xy = source_crds, 
                      mu = 0.08, 
                      md = 140, 
                      slp_thresh = 35, 
                      exp_div = 2.5, 
                      per_fct = 1.95, 
                      walks = 1000,
                      source_connect = TRUE,
                      connect_feature = feature_mask)

# Plot structure of sim_paths
str(sim_paths)
#> List of 4
#>  $ start_cell    : num 2803625
#>  $ cell_trav_freq: 'table' int [1:2027(1d)] 5 3 14 9 10 1 50 51 33 16 ...
#>   ..- attr(*, "dimnames")=List of 1
#>   .. ..$ all_indices: chr [1:2027] "2774166" "2775633" "2775634" "2775635" ...
#>  $ cell_max_vel  : num [1:2027] 15.2 12 13.9 16 17 ...
#>  $ prob_connect  : num 0.859

# Get connectivity probability
sim_paths$prob_connect
#> [1] 0.859

tp_raster <- walksToRaster(sim_paths, dem, method = "cdf_prob")
conn_raster <- connToRaster(sim_paths, dem)

# Crop results
tp_crop <- crop(tp_raster, ext(runout_polygon)+ 500)
conn_crop <- crop(conn_raster, ext(runout_polygon)+ 500)

# Plot results
par(mfrow = c(2,1))
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Traverse prob.")
plot(st_geometry(river_channel), col = "lightblue", add = TRUE)
plot(tp_crop, legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)


plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Connectivity prob.")
plot(st_geometry(river_channel), col = "lightblue", add = TRUE)
plot(conn_crop, legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)

Simulate runout from multiple-sources

To simulate runout from multiple source cell locations, we run runoutSim() iteratively, changing the xy coordinates during each iteration.

To create a list of xy coordinates from source points (vector), we use the makeSourceList() function.

# Get coordinates of source points to create a source list object
source_l <- makeSourceList(source_points)

# Perform random walk simulation for each source point
sim_runs <- list()

for(i in 1:length(source_l)){
  
  sim_runs[[i]] <- runoutSim(dem = dem, xy = source_l[[i]], 
                              mu = 0.08, 
                              md = 40, 
                              slp_thresh = 40, 
                              exp_div = 3, 
                              per_fct = 1.9, 
                              walks = 1000)
}

We can now convert this list of multiple random walks from different source cells to raster data using walksToRaster() and velocityToRaster().

# Coerce results to a raster
trav_freq <- walksToRaster(sim_runs, method = "freq", dem)
vel_ms <- velocityToRaster(sim_runs, dem, method = "max")
trav_prob <- walksToRaster(sim_runs, method = "max_cdf_prob", dem)


# Plot random walks from mulitple source cells
par(mfrow = c(1,3))
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), 
     main = "Traverse frequency")
plot(trav_freq, add = TRUE, alpha = 0.7)

plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), 
     main = "Traverse probability (max.)")
plot(trav_prob, add = TRUE, alpha = 0.7)

plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), 
     main = "Velocity (m/s)")
plot(vel_ms, col = map.pal("plasma"), add = TRUE, alpha = 0.7)

We can also view the results through an interactive leaflet map.


# Create interactive leaflet map
sim_map <- 
  # Debris flow observations
  leafmap(runout_polygons, opacity = 0.3, label = "Runout polygons") %>%
  leafmap(source_points, color = '#e74c3c', label = "Source points") %>%
  # Terrain data
  leafmap(hill, palette = grey(0:100/100), opacity = 1, add_legend = FALSE, 
            add_image_query = FALSE, label = "Hillshade") %>%
  leafmap(dem, palette = viridis::mako(10), opacity = 0.6, 
          add_image_query = FALSE, label = "DEM") %>%
  # Modelling results
  leafmap(trav_prob,  palette = viridis::viridis(10, direction = -1), 
          label = "Traverse prob", add_image_query = FALSE) %>%
  leafmap(vel_ms,  palette = viridis::plasma(10, direction = -1),
          label = "Velocity", add_image_query = FALSE)

# Start with mapped zoomed in
sim_map <- leaflet::setView(sim_map, lng = -70.13694, lat = -33.3703, zoom = 13)


# Hide layers
sim_map <- leaflet::hideGroup(sim_map, c("Velocity", "Source points", "DEM", "Hillshade"))

sim_map